## This code estimates the ideal points from the raw voting data
## This code takes a while to run on most machines as it uses MCMC
## Note that we provide the full posterior (Posterior.Rda) so it is possible to run diagnostics without doing the full estimation.


coalitions <- read.csv("coalitions.csv")

library(dplyr)
library(tidyr)
library(MCMCpack)
df <- tbl_df(coalitions)

#Filter out "nationals" and non-regular judges as well as judges who voted 10 times or less

df <- filter(df, !is.na(judgeCcode),
             judgeCcode != COWcode_1, 
             is.na(COWcode_2) |judgeCcode != COWcode_2,  
             is.na(COWcode_3) | judgeCcode != COWcode_3, 
             is.na(COWcode_4) | judgeCcode != COWcode_4, 
             is.na(COWcode_4) | judgeCcode != COWcode_5, 
             TotalJudgeVotes>10)

#Create Vote ID variable 
df$vid <- df %>% group_indices(filename, issue)
#create judge ID
df$jid <- df %>% group_indices(lastname)

RollCallList = unique(df$vid)	## CHECK: length(RollCallList)
JudgeList = unique(df$lastname)	## CHECK: length(JudgeList)

# (ii) Generate vote matrix
VoteMatrix = matrix(NA, nrow = length(JudgeList), ncol = length(RollCallList))
## dim(VoteMatrix)    ## Check dimensions of matrix	length(RollCallList)
colnames(VoteMatrix) = RollCallList
rownames(VoteMatrix) = JudgeList
## View(VoteMatrix)

R <-nrow(df)
# (iii) add values into the vote matrix
for (ii in 1: R){
  VoteMatrix[df$jid[ii],df$vid[ii]] <- df$progov[ii]
}


##Identification is to constrain all beta parameters to be positive

item.constraints.list  <- rep(list("colname" = list(2, "+")), times = ncol(VoteMatrix)) 
names(item.constraints.list) <- colnames(VoteMatrix)

posterior2   <- MCMCirtKdRob(VoteMatrix, dimensions=1,
                             burnin=200000, mcmc=1000000, thin=2000, verbose=50000,
                             B0=.25, b0=1, store.item=TRUE, store.ability=TRUE,
                             item.constraints = item.constraints.list,
                             ability.constraints=list("tulkens"=list(1,"-"),
                                                      "matscher"=list(1,"+")))

means.sds2 = summary(posterior2)[[1]][,1:2] # the first two columns are means and sds
theta2 = means.sds2[grepl("theta", rownames(means.sds2)),]

theta2 <- as.data.frame(theta2)
theta2$jid <- seq.int(nrow(theta2))

beta2 = as.data.frame(means.sds2[grepl("beta", rownames(means.sds2)),])


alpha2 = as.data.frame(means.sds2[grepl("alpha", rownames(means.sds2)),])


judges<-df %>%
  group_by(jid, lastname) %>%
  summarize(votes = n_distinct(vid)) 


judges <- left_join(judges, theta2, by = "jid")

load("JudgeInfo.RData")

judgeInfo <- full_join(judgeInfo, judges,  by = "lastname")

judgeInfo$theta<-judgeInfo$Mean



#Save ideal points
write.csv(judgeInfo, "Idealpointestimates.csv")
#Save posteriors

save(posterior2,file="posterior.Rda")
save(beta2,file="Discrimination parameters.RData")
save(alpha2,file="CutPoints.RData")


